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Abstract 

A method to derive the corrections for the dispersion of the reaction 
plane at intermediate energies is proposed. The method is based on 
the correlated, non-isotropic Gaussian approximation. It allowed to 
construct the excitation function of genuine flow values for the Au+Au 
reactions at 40-150 MeV/nucleon measured with the INDRA detector 
at GSI. 

1 Introduction 

Flow measurements are being performed since two decades [1, 2] starting 
from energies as low as few tens of MeV/nucleon up to the ultra relativistic 
energies available now, and for varying system sizes and asymmetries from 
C+C to Bi+Bi. Both, sideward collective motion (directed flow related ob- 
servables) and azimuthal anisotropics (elliptic flow related observables) have 
been extensively studied. The primary motivation to study these observables 
is their presumed link to the density and pressure attained during the heavy 
ion collisions. A precise knowledge of the behavior of the excitation function 
of both flow components is necessary to constrain the transport models in 
order to study the properties of nuclear matter far from the normal condi- 
tions. 

To cover a possibly large range of incident energies, different experiments 
using different detectors have to be merged. The results may depend on the 



acceptance and resolution of a given apparatus as well as on the method 
used to determine the flow observablcs and in particular the reaction plane. 
Thus merging different experiments requires representing their results in 
a manner least influenced by the instrumental and methodological factors. 
The instrumental factors are minimized by using high resolution and high 
coverage devices or correcting for inefficiencies. The method related ones, in 
case of the flow observables, are minimized by presenting the values corrected 
for the dispersion of the estimated reaction plane. 

Since the detectors can not measure the positions, spins and momenta of 
the reaction products simultaneously, the impact vector, and in particular 
the reaction plane, can only be estimated using the momenta, and the pre- 
cision of the flow measurement will depend on the accuracy of this estimate. 
Since the beginning of the flow investigation numerous methods have been 
proposed to determine the reaction plane dispersion [3-17], however all of 
them are based on the central limit theorem and thus limited to collisions 
at high energies which fulfill the high multiplicity requirement. 

The established methods proved their usefulness for correcting measured 
flow values at high energies (sec e.g. [3, 18-21]). They arc, however, not 
adequate for the intermediate energy reactions, especially below about 100 
AMeV. The present study was motivated to fill this gap. 



2 The method 

Originally, the directed flow has been quantified by measuring the in- 
plane component of the transverse momentum [3] and the elliptic flow by 
parametrizing the azimuthal asymmetries in terms of the "squeeze-out" ra- 
tio [22]. More recently, it has been proposed to express flow components 
in terms of the Fourier coefficients, vi, V2 and higher, obtained from the 
Fourier decomposition [9, 10, 12, 17] of the azimuthal distributions of the re- 
action products measured with respect to the reaction plane, with azimuthal 
angle (pR, and as a function of the particle type, rapidity y and possibly the 
transverse momentum p^- 
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The first two coefficients, v\ = (cos(0 — and = (cos2((/) — ^r)), 
characterize the directed and elliptic flow, respectively. 

Since the azimuth of the reaction plane, can only be estimated with a 
finite precision, the measured coefficients v^^"-^ are biased. They are related 



to the genuine ones through the following expression [10] : 

v~ = (cosn((/> - (/>£;)) = (cosn(</> -(pR + cpR- (^e)) = VnicosnAc/)) (2) 

where the average cosine of the azimuthal angle between the true and the 
estimated planes, (cosnAcj)) = {cos n{(j)R — (j)E)), is the required correction 
for a given harmonic. 

The reaction plane can be determined by several methods, including the 
flow-tensor method [23], the fission-fragment plane [24], the flow Q- vector 
method [3], the transverse momentum tensor [6] ("azimuthal correlation" 
[25]) method or others [26]. 

Taking advantage of the additivity of the Q-vector, which is defined 
as a weighted sum of the transverse momenta of the measured A'' reaction 
products: 

N 

Q = Y.o^iPi (3) 
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we adopt this method for estimating the reaction plane. The correction 
factor is then searched for using the sub-event method [3,10], which consists 
in splitting randomly each event into two equal multiplicity sub-events and 
getting the correction factor from the distribution of the relative azimuthal 
angle between the Q- vectors for sub-events ( "sub-Q- vectors" ) by fitting the 
theoretical distribution to it. 

For high multiplicity events measured at high energies this theoretical 
distribution can be obtained from the central limit theorem and assuming 
that the sub-Q-vectors are normally distributed. Ref. [10] gives an analytical 
formula for such a distribution for the case that the distributions of sub-Q- 
vectors arc independent and isotropic around their mean values. 

At low energies the particle multiplicities are lower and the events are 
characterized by a broad range of masses of the reaction products. Thus, 
the applicability of the central limit theorem for devising the corrections 
is less obvious. Nevertheless, taking advantage of the presumed relatively 
good memory of the entrance channel plane by the heavy remnants, it can 
be argued that the assumption of the Gaussian distribution of the flow Q- 
vector may still hold, at least in the first approximation. QMD-CHIMERA 
calculations, using a version of the code with a strict angular momentum 
conservation, confirm this assumption [27], and, even more, show that the 
random sub-events are also normally distributed even at 40 AMeV, except 
for very peripheral collisions. The last observation is crucial since the pro- 
posed method is based on the assumption of the normal distributions of 
Q-vectors for sub-events. The difference, as compared to the high energy 



case, is that the sub-events are strongly correlated [8] and that the distri- 
butions of the Q-vector are no longer isotropic but rather elongated in the 
reaction plane, as suggested in [10]. The elongation presumably reflects the 
increasing role of the in-plane emissions at low energies. 

Taking into account these observations, and following the method out- 
lined in appendix A of [17], we have derived the form of the joint probability 
distribution of the random sub-Q-vectors by imposing the constraint of mo- 
mentum conservation on the A^-particle transverse momentum distribution 
and by using the saddle-point approximation. 

The resulting distribution is a product of two bivariate Gaussians: 
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where we followed the convention of [10] of including the \/2 in a; the sub- 
scripts 1, 2 and s refer to the "sub-event". This distribution differs from 
those proposed in [8,10,17] in that it combines all three effects that influence 
the reaction plane dispersion at intermediate energies, namely the directed 
flow (through the mean in-plane component Qg or the resolution parameter 
Xs = Qs/o'sx [10]), the elliptic flow (through the ratio a = agx/o'sy) and the 
correlation between the sub-events [8] (through the correlation coefficient 

pe[-i,i]). 

Making the division into sub-events random ensures that the distribu- 
tions of the sub-Q-vectors are equivalent, in particular they have the same 
mean values and variances. Since the total-Q-vector is the sum of the sub- 
Q-vectors, Q = Qi + Q2, its distribution takes the following form (cf eq. 
(15) of ref. [10]): 
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From this, taking into account the definition of x [10] oiis finds that: 



X = Xsp/{l + p) (6) 

Note that, in principle, the sub events do not have to be of equal multi- 
plicity, however, keeping this condition reduces the fluctuations and makes 
the Gaussian assumption better fulfilled. 

As in [10] the use of the joint probability distribution Q is done by 
integrating it over the magnitudes of the Q-vectors and one angle, leav- 
ing the relative angle between the sub-events. Unlike in [10], the resulting 



distribution can not be presented in an analytical form, instead, it can be 
calculated numerically (two-dimensional integral) using e.g. the combined 
Gauss-Legendre/Gauss-Chebysliev quadrature. It depends on 3 parameters 
(Xs, a, p) which can be obtained from fitting it to the experimental or model 
data. 

The correction factors for the n-th harmonic fn, depending now on x and 
a, can be calculated (also numerically) as the mean values of the cosnA(/) 
obtained over the total-Q- vector distribution ©, in a similar way as in [10]. 




Figure 1: Correction factors for 
the first (3 upper curves) and 
the second (3 lower curves) har- 
monic as a function of the res- 
olution parameter, x, for differ- 
ent aspect ratios, a. 
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Figure n shows how the elongation of the Gaussian (a), or elliptic flow, 
modifies the correction factors for the first two harmonics. It demonstrates 
that the in-plane emissions (q > 1) enhance slightly the resolution for vi 
and considerably for V2 - even in the absence of the directed flow. On the 
other hand, squeeze-out (a < 1) deteriorates the resolution. In particular, it 
shows that the correction can be negative in case of small directed flow and 
squeeze-out. The relation ® between Xs and x indicates that the resolution 
improves in case the sub-events are anti-correlated (/> < 0). 

The introduction of 3 parameters in Q) was necessary to fit the whole 
family of shapes of the experimental distributions of the relative angle be- 
tween the sub-events. 

Fig. 121 shows an example of shapes that can be described by the inte- 
grated distribution (jlj. The parameters used for each curve are specified 
in the caption together with the corresponding correction factors. Case (a) 
corresponds to events with small directed flow and strong anti-correlation, 
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Figure 2: Distribution of rela- 
tive angle between the sub-Q- 
vectors for different parameters 
together with the resulting cor- 
rections. 
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producing a pronounced backward peaking observed at low energies and 
more central collisions. Case (b) corresponds to less central collisions at 
low energies. Case (e) corresponds to mid-central collisions at intermediate 
energies, where the directed flow increases. At higher energies the peak at 
small relative angles becomes narrower. 

Examples (c) and (d) deserve a special comment. They are almost indis- 
tinguishable. Case (d) represents an almost isotropic distribution with no 
directed flow and for independent sub-events (x = 0, p = 0) with the vanish- 
ing correction for directed flow ((cos Ac/)) = 0), while (c) represents the case 
with non-zero directed flow and moderate anti-correlation {p = —0.25). In 
this case the corrections are substantially different from zero. This example 
shows the difficulties that can be encountered while fitting the experimen- 
tal distributions. In these special cases the method may lose its resolving 
power, unless one can find a way to constrain some of the parameters, e.g. p. 
One way to do this is to get a hint on the value of the correlation coefficient 
from model calculations. For example, the CHIMERA calculations predict 
this coefficient to be around -0.43 for 40 AMeV and 2 < 6 < 8 fm and 
about -0.2 at 150 AMeV, indicating that the case (d) with vanishing flow, 
is less likely to represent the experimental distribution. Alternatively, one 
can impose constraints by calculating some rotational invariants from the 



measured data, e.g. {Qi • Q2) and (Qf) or (Q^)- From these the fohowing 
ratio can be constructed: 

{Qi ■ Q2) ^ . ^ p + a^ip + 2xl) , . 

{QD l + a2(l + 2x2) ^'^ 

This relation allows to reduce the number of fit parameters to two by 
using the experimentally accessible ratio f3. 

For the cases (c) and (d) from the previous example /? takes the values of 
-0.0114 and 0, respectively. Taking into account that f3 can be known with 
high accuracy for high statistics data (a few percent at worst), its precise 
knowledge constrains the fitting routine to search for a conditional minimum 
and helps to resolve the ambiguous cases. 

Instead of fitting the relative angle distributions one can express the 
probability distribution Q in terms of projections of one sub-Q- vector in the 
reference frame of the other. The corresponding 2-dimensional experimental 
distributions can then be fit using such a formula. The drawback might be 
that this formula depends on the magnitudes of the sub-Q- vectors, which 
introduces a dependence not only on the experimental uncertainties of the 
measured angles but also on the accuracy of the energy calibration. The 
advantage is that it can be calculated using only one-dimensional numerical 
integration. This formula contains four parameters, with agx being the 
additional one. However, using the calculated invariants, it is again possible 
to reduce the number of fit parameters to two. 



3 Preliminary results 

The experimental data unavoidably suffers from inaccuracies and lacking 
completeness. The fitting procedure yields relatively accurate results for 
the corrections for the first two harmonics in case of the simulations (e.g. 
2-5% for 40 AMeV and 0.2-0.4% for 150 AMeV and 4 < 6 < 8 fm), but 
in the case of the experimental data it will certainly return some effective 
correction factors biased by the experimental uncertainties. 

In order to minimize the experimental uncertainties one can either re- 
strict the data sample to "quasi-complete" events, with the possibility that it 
is no longer representative, or one can try to complete the event information 
by applying conservation laws, with the possibility of generating additional 
fluctuations. 

In the following, we adopt the latter approach. We analyze all the events 
which contain at least 35% of the total charge and complete them by adding 



a single "extra" fragment carrying the missing momentum, charge and mass. 
This procedure changes substantially the distributions of the relative angle 
between sub-events for peripheral collisions where, due to the energy thresh- 
olds, the heavy target-like fragment is always lost. The distributions become 
narrowly peaked at small relative angles which improves the resolution of 
the reaction plane. 
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Figure 3: Correction parameters for vi (left) and V2 (right) for 4 < & < 6 fm 
collisions of Au+Au measured with the INDRA detector. 

Fig. El presents the correction parameters for the first two harmonics 
calculated on the basis of the parameters Xs,Ci,P obtained from the fits to 
the experimental distributions of the relative angle between sub-events (solid 
circles), and from the fits of the projections of one sub-Q- vector on the other 
(solid squares). In both cases the calculated rotational invariants were used 
to constrain the fits. The stars represent the results of the fits with the 
fixed p as obtained from the CHIMERA calculation, and the open circles 
represent the results obtained using the standard Ollitrault method [10], 
which corresponds to the case of p = and a = 1. The figure demonstrates 
the utility and accuracy of the proposed method in cases where the standard 
methods do not apply, i.e. at low energies. It shows that the standard and 
the extended method approach each other at higher energies, but even at 
150 AMeV there is still about 5% and 15% discrepancy for the vi and V2 
corrections, respectively. Since at 150 AMeV the elliptic fiow is small, this 
discrepancy comes mostly from correlations between the sub-events which 
are still non-negligible (p ~ —0.2) at this energy. Differences between the fit 
results obtained using the angular distributions and the projections, based 
on the new method (solid symbols), refiect its systematic uncertainty which, 
apart from the 15 AMeV case, can be quantified as about 2% for directed 



flow and up to about 10% for the elliptic flow. 
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Figure 4: Corrected slopes of vi at midrapidity (left) and V2 at midrapidity (right) 
for Z—2 particles from mid-central (4 < < 6 fm) collisions of Au+Au measured 
whh the INDRA detector. 

Figure shows the excitation function of the slope of vi at midrapidity 
(left) and of V2 at midrapidity (right) corrected for the dispersion of the 
reconstructed reaction plane using the corrections obtained with the present 
method (angular distribution fits). 

The vi and V2 parameters have been calculated in the beam frame, ex- 
cluding the particle of interest from the reaction plane reconstruction and 
correcting for the momentum conservation [17]. The linear fits to obtain the 
slope parameters have been carried out in the range of ±0.4 of the scaled 
CM rapidity for all energies except 15 AMeV, where the range of ±0.55 has 
been used to improve the accuracy. In the case of the second harmonic the 
value of V2 at midrapidity has been obtained by fitting a parabola to the V2 
vs scaled rapidity relation and taking the value of the fit at ycM = 0. The 
point at 15 AMeV has been included in the systematics to get a hint on the 
evolution of the excitation function at very low energies. Due to the very 
low statistics, this result should be taken with care and, possibly, verified 
with dedicated and more complete measurements. 

The trends observed for the uncorrected data [28] for vi are preserved. 
The excitation function for Z=2 changes sign between 50 and 60 AMeV and 
continues to decrease at lower energies. 

The excitation function of the elliptic fiow shows a possible saturation 
of the in-plane enhancement near 0.08 for bombarding energies below about 
60 AMeV or even a signature of a maximum near 50 AMeV. 



4 Summary 



We have proposed a new method to calculate the corrections resulting from 
the dispersion of the reaction plane at intermediate energies. The method, 
which in some respect is an extension of the one proposed by Ollitrault [10], 
takes into account the combined effects of the directed flow, elliptic flow 
and correlations on the reaction plane resolution. It shows an increasing 
role of correlations and of the in-plane enhancement at lower energies. The 
corrected experimental excitation functions of vi and V2 constitute a pre- 
cise constraint for the dynamical models which are aimed at inferring the 
properties of nuclear matter at varying densities and pressures from flow 
observables. 
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